We share this use case to show all features of rANOMALY package that allow to highlight microbial community differences between groups of samples. rANOMALY main features are :

  • ASV (Amplpicon Sequence Variant) computation thanks to dada2 algorithm

  • Taxonomic assignment thanks to IDTAXA from DECIPHER package. rANOMALY allow to assign with up to 2 complementary reference databases (eg. SILVA which is generalist, DAIRYdb which is specific to dairy environments), keep the best and more confident taxonomic assignment, check for taxonomy incongruencies like empty fields or multiple ancestors.

  • phyloseq format to handle data and downstream statistical analysis.

  • Decontamination step thanks to decontam package. rANOMALY allows to filter ASV depending on ASV overall frequency (low abundance filter), and ASV prevalence (presence in 1 or more samples), specific filters like manual suppression of specific taxa are available to.

  • Complete statistical analysis workflow: rarefaction curves, community composition plots, alpha/beta diversity analysis with associated statistical tests and multivariate analysis ; differential analysis using four different methods (metacoder, DESeq2, metagenomeSeq, plsda).

  • At each step, plots and R objects are returned in command line AND saved in folders to be able to relaunch the analysis and to export plots.

1 Help

Each function has a detailed help accessible in R via ?{function}.

2 Tests datasets

The dataset can be downloaded via this link.

This tutorial assumes that you have extracted all the read file in a folder named reads along with the sample-metadata.csv file.

We share a 24 DNA samples test dataset extracted from rats feces at two different time (t0 & t50) and in two nutrition conditions. Rats were fed with two strains of some bacteria (wild type & mutant). Also, extraction control sample (blank) are included.

sm <- read.table("sample_metadata.csv", sep="\t",header=TRUE)
DT::datatable(sm)
load("decontam_out/robjects.Rdata")

3 Processing of raw sequences

3.1 ASV definition with DADA2

The first step is the creation of ASVs thanks to the dada2 package. In rANOMALY, only one function is needed to compute all the different steps require from this package. The function is designed to process Illumina paired end sequences in fastq format.

Sample names are extracted from the file name, from the beginning to the first underscore (_), so files must be formatted as followed : {sample-id1}_R1.fastq.gz {sample-id1}_R2.fastq.gz etc… Those must match the sample names stored in sample-metadata.csv file.

dada_res = dada2_fun(path="./reads", dadapool = "pseudo", compress=TRUE, plot=FALSE)

Main outputs:

  • ./dada2_out/read_tracking.csv summarizes the read number after each filtering step.

  • ./dada2_out/raw_otu-table.csv the raw ASV table.

  • ./dada2_out/rep-seqs.fna: fasta file with all representative sequences for each ASV.

  • ./dada2_out/robjects.Rdata with saved dada_res list containing raw ASV table and representative sequences in objects otu.table, seqtab.export & seqtab.nochim.

Read tracking table:

  • input: raw read number.

  • filtered: after dada2 filtering step: no N’s in sequence, low quality, and phiX.

  • denoisedF & denoisedR: after denoising. Forward & Reverse.

  • merged: after merging R1 & R2.

  • nonchim: after chimeras filtering.

DT::datatable(read.table("dada2_out/read_tracking.csv",sep="\t",header=TRUE), filter = "top")

3.2 Taxonomic assignment

assign_taxo_fun function uses IDTAXA function from DECIPHER package, and allows to use 2 different databases. It keeps the best assignment on 2 criteria, resolution (depth in taxonomy assignment) and confidence. The final taxonomy is validated for multiple ancestors and incongruities correction step.

We share the latest databases we use in the IDTAXA format in this link. You can also generate your own IDTAXA formatted database following those instructions and scripts we provide at this page.

tax.table = assign_taxo_fun(dada_res = dada_res, id_db = c("path_to_your_banks/silva/SILVA_SSU_r132_March2018.RData","path_to_your_banks/DAIRYdb_v1.2.0_20190222_IDTAXA.RData") )

Main file outputs:

  • ./idtaxa/robjects.Rdata with taxonomy in phyloseq format in tax.table object.

  • final_tax_table.csv the final assignation table that will be use in next steps.

  • allDB_tax_table.csv raw assignations from the two databases, mainly for debugging.

3.3 Phylogenetic Tree

The phylogenetic tree from the representative sequences is generated using phangorn and DECIPHER packages.

tree = generate_tree_fun(dada_res)

Main outputs: - tree_robjects.Rdata with phylogenetic tree object in phyloseq format.

3.4 Phyloseq object

To create a phyloseq object, we need to merge four objects and one file:

  • the asv table dada_res and the representative sequences seqtab.nochim from dada2_robjects.Rdata

  • a taxonomy table (taxtable) taxo_robjects.Rdata from taxo_robjects.Rdata

  • the phylogenetic tree tree from tree_robjects.Rdata

  • metadata from sample-metadata.csv

data = generate_phyloseq_fun(dada_res = dada_res, tax.table = tax.table, tree = tree, metadata = "./sample_metadata.csv")
data
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 146 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 5 sample variables ]
## tax_table()   Taxonomy Table:    [ 146 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 146 tips and 144 internal nodes ]
## refseq()      DNAStringSet:      [ 146 reference sequences ]

Main output:

  • ./phyloseq/robjects.Rdata with phyloseq object in data for raw counts and data_rel for relative abundance.

3.5 Decontamination

The decontam_fun function uses decontam R package with control samples to filter contaminants. The decontam package offers two main methods, frequency and prevalence (and then you can combine those methods). For frequency method, it is mandatory to have the DNA concentration of each sample in phyloseq (and hence in the sample-metadata.csv). The prevalence method does not need DNA quantification, this method allows to compare presence/absence of ASV between real samples and control samples and then identify contaminants.

Tips: sequencing plateforms often quantify the DNA before sequencing, but do not automaticaly give the information. Just ask for it ;).

Our function integrates the basic ASV frequency (nb_reads_ASV/nb_total_reads) and overall prevalence (nb_sample_ASV/nb_total_sample) filtering. We included an option to filter out ASV based on their taxa names (for known recurrent contaminants).

data = decontam_fun(data = data, domain = "Bacteria", column = "type", ctrl_identifier = "control", spl_identifier = "sample", number = 100, krona= TRUE)

Main outputs: - robjects.Rdata with contaminant filtered phyloseq object named data. - Exclu_out.csv list of filtered ASVs for each filtering step. - Krona plot before and after filtering. - raw_asv-table.csv & relative_asv-table.csv. - venndiag_filtering.png.

4 Plots, diversity and statistics

We are currently developping a ShinyApp to visualize your data, sub-select your samples/taxons and do all those analyses interactively !!! ExploreMetabar

4.1 Rarefaction curves

In order to observe the sampling depth of each samples we start by plotting rarefactions curves. Those plots are generated by plotly which makes the plot interactive.

rareplot = rarefaction(data, "strain_time", 100 )
htmltools::tagList(list(rareplot))

4.2 Composition plots

The bar_fun function outputs a composition plot, in this example it reveals the top 10 genus present in our samples. The function allows to plot at different rank and to modify the number of taxas to show.

  • Ord1 option order the sample along the X axis.

  • Fact1 option control labels of the X axis. Fact1="sample.id" if you don’t want the sample to be renamed.

4.2.1 Raw abundance

bars_fun(data = data, top = 10, Ord1 = "strain_time", Fact1 = "strain_time", rank="Genus", relative = FALSE)

4.2.2 Relative abundance

bars_fun(data = data, top = 10, Ord1 = "strain_time", Fact1 = "strain_time", rank="Genus", relative = TRUE)

4.3 Diversity analyses

4.3.1 Alpha diversity

This function computes various alpha diversity indexes, it uses the estimate_richness function from phyloseq.

Available measures : c(“Observed”, “Chao1”, “ACE”, “Shannon”, “Simpson”, “InvSimpson”, “Fisher”).

It returns a list which contains:

  • a boxplot comparing conditions. ($plot)

  • a table of indices values. ($alphatable)

And for each of the computed indices :

  • an ANOVA analysis. (${measure}$anova)

  • a pairwise wilcox test result comparing conditions and giving the pvalue of each comparison tested. (${measure}$wilcox_col1, ${measure}$wilcox_col2_fdr, ${measure}$wilcox_col2_collapsed)

  • a mixture model if your data include repetition in sampling, ie. column3 option. (${measure}$anovarepeat, ${measure}$mixedeffect)

All this in a single function.

alpha <- diversity_alpha_fun(data = data, output = "./plot_div_alpha/", column1 = "strain", column2 = "time", column3 = "", supcovs = "", measures = c("Observed", "Shannon") )

4.3.2 Table of diversity index values

The table of values for each indices you choose to compute.

pander(alpha$alphatable, style='rmarkdown')
  Observed Shannon
SB1.Sauv0 41 1.477
SB10.Mut0 40 2.073
SB11.Mut0 51 2.178
SB12.Mut0 38 2.116
SB13.Sauv50 46 2.691
SB14.Sauv50 57 2.905
SB15.Sauv50 50 2.793
SB16.Sauv50 52 2.8
SB17.Sauv50 49 2.624
SB18.Sauv50 54 2.831
SB19.Mut50 66 2.638
SB2.Sauv0 26 2.099
SB20.Mut50 72 2.721
SB21.Mut50 79 3.062
SB22.Mut50 81 2.81
SB23.Mut50 84 3.175
SB24.Mut50 90 3.148
SB3.Sauv0 19 0.1962
SB4.Sauv0 41 2.52
SB5.Sauv0 46 1.923
SB6.Sauv0 46 1.067
SB7.Mut0 33 2.256
SB8.Mut0 58 2.089
SB9.Mut0 50 2.237

4.3.3 Boxplots

The boxplots of those values.

ggplotly(alpha$plot) %>%
  layout(boxmode = "group")
## Warning: 'layout' objects don't have these attributes: 'boxmode'
## Valid attributes include:
## 'font', 'title', 'uniformtext', 'autosize', 'width', 'height', 'margin', 'paper_bgcolor', 'plot_bgcolor', 'separators', 'hidesources', 'showlegend', 'colorway', 'datarevision', 'uirevision', 'editrevision', 'selectionrevision', 'template', 'modebar', 'meta', 'transition', '_deprecated', 'clickmode', 'dragmode', 'hovermode', 'hoverdistance', 'spikedistance', 'hoverlabel', 'selectdirection', 'grid', 'calendar', 'xaxis', 'yaxis', 'ternary', 'scene', 'geo', 'mapbox', 'polar', 'radialaxis', 'angularaxis', 'direction', 'orientation', 'editType', 'legend', 'annotations', 'shapes', 'images', 'updatemenus', 'sliders', 'colorscale', 'coloraxis', 'metasrc', 'barmode', 'bargap', 'mapType'

4.3.4 Tests on Observed index

4.3.4.1 ANOVA results

For each indices, you have access to the ANOVA test. Here we present the result for the “Observed” indice.

pander(alpha$Observed$anova)
Analysis of Variance Model
  Df Sum Sq Mean Sq F value Pr(>F)
Depth 1 49.36 49.36 0.5091 0.4838
strain 1 1877 1877 19.36 0.0002764
time 1 3649 3649 37.64 5.392e-06
Residuals 20 1939 96.96 NA NA

4.3.4.2 Wilcox test

Wilcox tests are made on each factor you have entered, and the combination of the two. Here “strain” and “time”.

4.3.4.3 Wilcox test for “strain” factor

pander(alpha$Observed$wilcox_col1)
  mutant
wildtype 0.043

4.3.4.4 Wilcox test for “time” factor

pander(alpha$Observed$wilcox_col2_fdr)
  t0
t50 0.001

4.3.4.5 Wilcox test for the collapsed factors

pander(alpha$Observed$wilcox_col2_collapsed)
  mutant_t0 mutant_t50 wildtype_t0
mutant_t50 0.002 NA NA
wildtype_t0 0.377 0.005 NA
wildtype_t50 0.336 0.002 0.008

4.4 Beta diversity

The diversity_beta_fun function is an exhaustive function allowing to generate all possible tests and ordination in one command.

The diversity_beta_light function is equivalent to the first function but allows to generate specific tests and figures ready to publish in rmarkdown as in the example below. It is based on the vegan package function vegdist for the distance calculation and ordinate for the ordination plot.

We include statistical tests to ease the interpretation of your results. A permutational ANOVA to compare groups and test that centroids and dispersion of the groups are equivalent for all groups. User can inform col and cov arguments to assess PERMANOVA to determine significant differences between groups (eg. factor “strain” here and covariable “time”). A pairwise-PERMANOVA is processed to determine which condition is significantly different from another (based on p-value).

As return, you will get an object that contains:

  • An ordination plot ($plot)

  • The permANOVA results ($permanova)

  • The pairwise permANOVA ($pairwisepermanova)

# beta <- diversity_beta_light(data = data, output = "./plot_div_beta/", glom = "ASV", column1 = "time", column2 = "strain", covar ="")
beta_strain = diversity_beta_light(data, col = "strain", cov="time", dist0 = "bray", ord0 = "MDS", output="./plot_div_beta_strain/", tests = TRUE)
beta_time = diversity_beta_light(data, col = "time", cov="strain", dist0 = "bray", ord0 = "MDS", output="./plot_div_beta_time/", tests = TRUE)

beta_strain_time = diversity_beta_light(data, col = "strain_time", dist0 = "bray", ord0 = "MDS", output="./plot_div_beta_strain_time/", tests = TRUE)

4.4.1 Ordination plot

htmltools::tagList(list(ggplotly(beta_strain$plot), ggplotly(beta_time$plot), ggplotly(beta_strain_time$plot)))

4.4.2 Statistical tests

  • permanova
pander(beta_strain$permanova)
Permutation: free
  Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
Depth 1 0.5075 0.5075 3.122 0.06845 0.01399
time 1 2.185 2.185 13.44 0.2946 0.000999
strain 1 1.471 1.471 9.049 0.1984 0.000999
Residuals 20 3.251 0.1626 NA 0.4385 NA
Total 23 7.415 NA NA 1 NA
  • pairwise permanova on concatenated factors
pander(beta_strain_time$pairwisepermanova)
Table continues below
pairs Df SumsOfSqs F.Model R2 p.value
wildtype_t0 vs mutant_t0 1 0.9528 4.64 0.317 0.007
wildtype_t0 vs wildtype_t50 1 2.021 28.97 0.7434 0.004
wildtype_t0 vs mutant_t50 1 2.197 26 0.7223 0.003
mutant_t0 vs wildtype_t50 1 1.681 11.59 0.5369 0.004
mutant_t0 vs mutant_t50 1 1.57 9.826 0.4956 0.003
wildtype_t50 vs mutant_t50 1 1.818 75.22 0.8827 0.001
p.adjusted sig
0.007 *
0.0048 *
0.0048 *
0.0048 *
0.0048 *
0.0048 *

4.5 Differential analyses

We choose three different methods to process differential analysis which is a key step of the workflow. The main advantage of the use of multiple methods is to cross validate differentially abundant taxa between tested conditions.

4.5.1 Metacoder

Metacoder is the most simple differential analysis tool of the three. Counts are normalized by total sum scaling to minimize the sample sequencing depth effect and it uses a Kruskal-Wallis test to determine significant differences between sample groups. The metacoder_fun function allows the user to choose the taxonomic rank, which factor to the test (column1), and a specific pairwise comparison (comp) to launch the differential analysis.

It also produces pretty graphical trees, representing taxas present in both groups and coloring branches depending on the group in which this taxa is more abundant. Two of those trees are produced, the raw one and the non significant features filtered one (p-value <= 0.05).

out1 = metacoder_fun(data = data, output = "./metacoder", column1 = "strain_time", rank = "Family", signif = TRUE, plottrees = TRUE, min ="10", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
  • Table
 DT::datatable(out1$table, filter = "top", options = list(scrollX = TRUE))
  • Comparaison 1
 out1$wildtype_t0_vs_mutant_t0$signif
## [[1]]

  • Comparaison 2
 out1$wildtype_t50_vs_mutant_t50$signif
## [[1]]

4.5.2 DESeq2

DESeq2 is a widely used method, primarily for RNAseq applications, for assessing differentially expressed genes between controlled conditions. It use for metabarcoding datas is sensibly the same. The deseq2_fun allows to process differential analysis as metacoder_fun, and users can choose the taxonomic rank, factor to test and which conditions to compare. DESeq2 algorithm uses negative binomial generalized linear models with VST normalization (Variance Stabilizing Transformation).

Main output is a list with :

  • list\({comparison}\)plot: plot with significant features
  • list\({comparison}\)table: table with statistics (LogFoldChange, pvalue, adjusted pvalue…)
#fig.keep='all', fig.align='left', fig.width = 15, fig.height = 10
out2 = deseq2_fun(data = data, output = "./deseq/", column1 = "strain_time", verbose = 1, rank = "Family", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(out2$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(out2$wildtype_t50_vs_mutant_t50$table, filter = "top", options = list(scrollX = TRUE))

4.5.3 MetagenomeSeq

MetagenomeSeq uses a normalization method able to control for biases in measurements across taxonomic features and a mixture model that implements a zero-inflated Gaussian distribution to account for varying depths of coverage. As deseq2_fun, metagenomeseq_fun returns a table with statistics and a plot with significant features for each comparison.

out3 = metagenomeseq_fun(data = data, output = "./metagenomeseq/", column1 = "strain_time", verbose = 1, rank = "Family", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(out3$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(out3$wildtype_t50_vs_mutant_t50$table, filter = "top", options = list(scrollX = TRUE))

4.5.4 Aggregate methods

The aggregate_fun function allows to merge the results from the three differential analysis methods used before, to obtain one unique table with all informations of significant differentially abundant features.

The generated table include the following informations :

  • seqid: ASV ID

  • Comparaison: Tested comparison

  • Deseq: differentially abundant with this method (0 no or 1 yes)

  • metagenomeSeq: differentially abundant with this method (0 no or 1 yes)

  • metacoder: differentialy abundant with this method (0 no or 1 yes)

  • sumMethods: sum of methods in which feature is significant.

  • DESeqLFC: Log Fold Change value as calculated in DESeq2

  • absDESeqLFC: absolute value of Log Fold Change value as calculated in DESeq2

  • MeanRelAbcond1: Mean relative abundance in condition 1

  • MeanRelAbcond2: Mean relative abundance in condition 2

  • Condition: in which the mean feature abundance is higher.

  • Taxonomy & representative sequence.

resF = aggregate_fun(data = data, metacoder = "./metacoder/metacoder_signif_Family.csv", deseq = "./deseq/", mgseq = "./metagenomeseq/", output = "./aggregate_diff/", column1 = "strain_time", column2 = NULL, verbose = 1, rank = "Genus", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(resF$wildtype_t0_vs_mutant_t0$plot)
ggplotly(resF$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(resF$table, filter = "top", options = list(scrollX = TRUE))

4.5.5 PLSDA (Partial Least Squares Discriminant Analysis)

Our PLS-DA/sparse PLS-DA analysis is based on the mixOmics package. These are supervised classification methods. They allow to define most important features discriminating our groups. The plsda_res function outputs a list of graphs and tables:

  • plsda$plotIndiv ordination plot from the PLSDA analysis.

  • splsda$plotIndiv ordination plot from the sparse PLSDA analysis.

  • loadings$comp{n} all the loadings plots for each component from the sPLSDA analysis. Loadings weights are displayed and colors correspond to the group in which each feature has its maximum mean value.

  • splda.loading_table Loading of all taxas for each component.

  • splsda$plotArrow arrow plot of samples.

plsda_res <- plsda_fun(data = data, output = "./plsda_family/", column1 = "strain_time", rank = "Family")
plsda_res$splsda.plotIndiv$graph

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp1.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp2.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp3.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp4.png')

5 Miscellaneous function

5.1 Heatmap

heatmap_fun allows to plot relative abundance of the top abundant taxas (20 by default). Users can choose which taxonomic rank to plot and the factor used to separate plots.

heatmap_plot = heatmap_fun(data = data, column1 = "strain_time", top = 20, output = "./plot_heatmap/", rank = "Species")
ggplotly(heatmap_plot)

5.2 Shared taxa

  • Venn Diagram

ASVenn_fun allows to define shared taxa between specific conditions (up to 5 conditions). It generates a venn diagram and a table informing presence of each ASV in each condition with its taxonomy and sequence. Counts can be verified with this table. rank argument allows to generate output with desired taxonomic level.

outvenn = ASVenn_fun(data = data,output = "./ASVenn/", rank = "ASV", column1 = "strain_time", shared = TRUE)
  • Venn diagram:
grid.draw(outvenn$venn_plot)

  • Table:
DT::datatable(outvenn$TABf, filter = "top", options = list(scrollX = TRUE))
  • Cytoscape phy2cyto_fun allows to generate input files (SIF format) for Cytoscape which is useful to visualize shared taxa. You can see below an example of representation with Cytoscape on our test data. Each nodes and arrows position and aesthetic can be easily modified.

5.3 Supplementary tools and features

csv2phyloseq_fun allows to import the 4 files (ASV, taxonomy, tree and sequences) in tabulated format to generate phyloseq object ready to use with rANOMALY.

assign_fasta_fun is used to assign sequences from any fasta files.

export_to_stamp_fun allows you to generate input files for STAMP.

split_table_fun allows to split the phyloseq object in multiple sub-object according to one factor (column).

update_metadata_fun allows you to easily modify sample_data part of the phyloseq object.

6 R Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ranomaly_0.0.0.9000 pander_0.6.3        DT_0.15            
## 
## loaded via a namespace (and not attached):
##   [1] taxa_0.3.4                  tidyselect_1.1.0           
##   [3] RSQLite_2.2.0               AnnotationDbi_1.48.0       
##   [5] htmlwidgets_1.5.1           grid_3.6.3                 
##   [7] BiocParallel_1.20.1         Rtsne_0.15                 
##   [9] devtools_2.3.2              IHW_1.14.0                 
##  [11] munsell_0.5.0               codetools_0.2-16           
##  [13] withr_2.2.0                 colorspace_1.4-1           
##  [15] Biobase_2.46.0              phyloseq_1.30.0            
##  [17] knitr_1.29                  rstudioapi_0.11            
##  [19] stats4_3.6.3                ggsignif_0.6.0             
##  [21] labeling_0.3                slam_0.1-47                
##  [23] GenomeInfoDbData_1.2.2      lpsymphony_1.14.0          
##  [25] hwriter_1.3.2               farver_2.0.3               
##  [27] bit64_0.9-7.1               rhdf5_2.30.1               
##  [29] rprojroot_1.3-2             vctrs_0.3.2                
##  [31] generics_0.0.2              lambda.r_1.2.4             
##  [33] xfun_0.16                   R6_2.4.1                   
##  [35] GenomeInfoDb_1.22.1         GA_3.2                     
##  [37] locfit_1.5-9.4              bitops_1.0-6               
##  [39] microbiome_1.8.0            DelayedArray_0.12.3        
##  [41] assertthat_0.2.1            scales_1.1.1               
##  [43] gtable_0.3.0                phangorn_2.5.5             
##  [45] processx_3.4.3              rlang_0.4.7                
##  [47] genefilter_1.68.0           splines_3.6.3              
##  [49] lazyeval_0.2.2              rstatix_0.6.0              
##  [51] broom_0.7.0                 yaml_2.2.1                 
##  [53] reshape2_1.4.4              abind_1.4-5                
##  [55] crosstalk_1.1.0.1           backports_1.1.8            
##  [57] tools_3.6.3                 usethis_1.6.3              
##  [59] ggplot2_3.3.2               ellipsis_0.3.1             
##  [61] gplots_3.0.4                decontam_1.6.0             
##  [63] biomformat_1.14.0           RColorBrewer_1.1-2         
##  [65] BiocGenerics_0.32.0         sessioninfo_1.1.1          
##  [67] Rcpp_1.0.5                  plyr_1.8.6                 
##  [69] zlibbioc_1.32.0             psadd_0.1.2                
##  [71] purrr_0.3.4                 RCurl_1.98-1.2             
##  [73] ps_1.3.3                    prettyunits_1.1.1          
##  [75] ggpubr_0.4.0                Wrench_1.4.0               
##  [77] viridis_0.5.1               cowplot_1.0.0              
##  [79] S4Vectors_0.24.4            SummarizedExperiment_1.16.1
##  [81] haven_2.3.1                 cluster_2.1.0              
##  [83] fs_1.4.2                    DECIPHER_2.14.0            
##  [85] magrittr_1.5                RSpectra_0.16-0            
##  [87] data.table_1.13.0           futile.options_1.0.1       
##  [89] pairwiseAdonis_0.0.1        openxlsx_4.1.5             
##  [91] ranacapa_0.1.0              matrixStats_0.56.0         
##  [93] pkgload_1.1.0               hms_0.5.3                  
##  [95] evaluate_0.14               xtable_1.8-4               
##  [97] XML_4.0-0                   VennDiagram_1.6.20         
##  [99] rio_0.5.16                  jpeg_0.1-8.1               
## [101] readxl_1.3.1                gridExtra_2.3              
## [103] IRanges_2.20.2              shape_1.4.4                
## [105] testthat_2.3.2              compiler_3.6.3             
## [107] ellipse_0.4.2               tibble_3.0.3               
## [109] KernSmooth_2.23-17          crayon_1.3.4               
## [111] htmltools_0.5.0             venn_1.9                   
## [113] corpcor_1.6.9               mgcv_1.8-33                
## [115] tidyr_1.1.0                 geneplotter_1.64.0         
## [117] RcppParallel_5.0.2          DBI_1.1.0                  
## [119] formatR_1.7                 MASS_7.3-53                
## [121] ShortRead_1.44.3            Matrix_1.2-18              
## [123] ade4_1.7-15                 car_3.0-8                  
## [125] permute_0.9-5               cli_2.0.2                  
## [127] quadprog_1.5-8              gdata_2.18.0               
## [129] parallel_3.6.3              igraph_1.2.5               
## [131] GenomicRanges_1.38.0        forcats_0.5.0              
## [133] pkgconfig_2.0.3             GenomicAlignments_1.22.1   
## [135] foreign_0.8-76              plotly_4.9.2.1             
## [137] foreach_1.5.0               rARPACK_0.11-0             
## [139] annotate_1.64.0             admisc_0.9                 
## [141] multtest_2.42.0             XVector_0.26.0             
## [143] stringr_1.4.0               callr_3.4.4                
## [145] digest_0.6.25               vegan_2.5-6                
## [147] dada2_1.14.1                Biostrings_2.54.0          
## [149] fastmatch_1.1-0             rmarkdown_2.3              
## [151] cellranger_1.1.0            curl_4.3                   
## [153] Rsamtools_2.2.3             gtools_3.8.2               
## [155] lifecycle_0.2.0             nlme_3.1-149               
## [157] jsonlite_1.7.1              Rhdf5lib_1.8.0             
## [159] mixOmics_6.10.9             carData_3.0-4              
## [161] futile.logger_1.4.3         viridisLite_0.3.0          
## [163] desc_1.2.0                  limma_3.42.2               
## [165] fansi_0.4.1                 pillar_1.4.6               
## [167] metacoder_0.3.4             lattice_0.20-41            
## [169] httr_1.4.2                  plotrix_3.7-8              
## [171] pkgbuild_1.1.0              survival_3.1-12            
## [173] glue_1.4.1                  remotes_2.2.0              
## [175] zip_2.0.4                   fdrtool_1.2.15             
## [177] png_0.1-7                   iterators_1.0.12           
## [179] glmnet_4.0-2                bit_1.1-15.2               
## [181] stringi_1.4.6               metagenomeSeq_1.28.2       
## [183] ggfittext_0.9.0             blob_1.2.1                 
## [185] DESeq2_1.29.13              latticeExtra_0.6-29        
## [187] caTools_1.18.0              memoise_1.1.0              
## [189] dplyr_1.0.0                 ape_5.4